function [datem,pcorrect,pcorrectRaw] = InverseBarometerCorrection(metfolder,yeari,NDBCname)
% InverseBarometerCorrection.m
% SNG 24 May 2012
% uses meterological pressure data to determine the proper inverse
% barometer correction for tidal gauges
% met data is from NBDC buoys: http://www.ndbc.noaa.gov/
% 
% SNG 18 Nov 2012 updated to be general for any year 2000 to 2010, to add
% more years, just download the .txt file at the above website for both
% Tatoosh Island (ttiw1) and Newport (nwpo3). Note there are other
% intermediate gauges that one could also include to be even more accurate


if nargin == 1
    NDBCname = 'ttiw1';
end

%choose the file based on the station and year
%tatoosh island
if strcmp(NDBCname,'ttiw1')
    metfile = ['ttiw1h' num2str(yeari) '.txt'];
%newport
elseif strcmp(NDBCname,'nwpo3')
    metfile = ['nwpo3h' num2str(yeari) '.txt'];
end
metin = [metfolder metfile];

%open the file
fid = fopen(metin);

%read in the header info
c2 = textscan(fid,'%s',20);

%read in the data to a cell array
%first re-open to read from beginning
fid = fopen(metin);
%do not read in the header line
extrahdr = 1;
disp(['reading ', metfile]); disp(' ');
c = textscan(fid,'%n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n','headerlines',extrahdr);
%close the raw data file
fclose(fid);

%extract times and pressure fields
year = c{1};
mo = c{2};
dy = c{3};
hr = c{4};
%find the time stamp of the atmospheric pressure field
bari = find(strcmp(c2{1},'BAR')==1);
if isempty(bari) % R. McCabe 07/2014 added since .txt files change format in 2007 (need to test!) 
     bari = find(strcmp(c2{1},'PRES')==1);
end
phPa = c{bari};

%date
datem = datenum([year mo dy hr zeros(size(hr)) zeros(size(hr))]); % R. McCabe buyer beware... not always the case
%year day
%jd = datem-datem(1); % R. McCabe commented out 07/2014 to avoid jd.
% R. McCabe 07/2014 also update the function call to have 'datem' instead of 'jd'
% [datem,pcorrect,pcorrectRaw] = InverseBarometerCorrection

%remove bad data
phPa(phPa>9000) = nan;

%fill in missing data with a linear interpolation
k = isnan(phPa);
%phPai = interp1(jd(~k),phPa(~k),jd); % R. McCabe commented out 07/2014 to avoid jd
phPai = interp1(datem(~k),phPa(~k),datem);

%godin low pass filter
%phPalp = Z_dasfilt(phPai,'godin');
phPalp = Z_godin(phPai); % R. McCabe 07/2014, Z_dasfilt is inferior junk!

%remove mean & convert from cm to m
pcorrectRaw = (phPai - nanmean(phPai))./100;
pcorrect = (phPalp - nanmean(phPalp))./100;



